home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / torture < prev    next >
Text File  |  1994-08-04  |  28KB  |  1,058 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26. /*
  27.     This file contains a series of tests for the Meschach matrix
  28.     library, parts 1 and 2
  29. */
  30.  
  31. static char rcsid[] = "$Id: torture.c,v 1.5 1994/05/18 01:53:58 des Exp $";
  32.  
  33. #include    <stdio.h>
  34. #include    <math.h>
  35. #include    "matrix2.h"
  36. #include        "matlab.h"
  37.  
  38. #define    errmesg(mesg)    printf("Error: %s error: line %d\n",mesg,__LINE__)
  39. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  40.  
  41. static char *test_err_list[] = {
  42.    "unknown error",            /* 0 */
  43.    "testing error messages",        /* 1 */
  44.    "unexpected end-of-file"        /* 2 */
  45. };
  46.  
  47.  
  48. #define MAX_TEST_ERR   (sizeof(test_err_list)/sizeof(char *))
  49.  
  50. /* extern    int    malloc_chain_check(); */
  51. /* #define MEMCHK() if ( malloc_chain_check(0) ) \
  52. { printf("Error in malloc chain: \"%s\", line %d\n", \
  53.      __FILE__, __LINE__); exit(0); } */
  54. #define    MEMCHK() 
  55.  
  56. /* cmp_perm -- returns 1 if pi1 == pi2, 0 otherwise */
  57. int    cmp_perm(pi1, pi2)
  58. PERM    *pi1, *pi2;
  59. {
  60.     int        i;
  61.  
  62.     if ( ! pi1 || ! pi2 )
  63.     error(E_NULL,"cmp_perm");
  64.     if ( pi1->size != pi2->size )
  65.     return 0;
  66.     for ( i = 0; i < pi1->size; i++ )
  67.     if ( pi1->pe[i] != pi2->pe[i] )
  68.         return 0;
  69.     return 1;
  70. }
  71.  
  72. /* px_rand -- generates sort-of random permutation */
  73. PERM    *px_rand(pi)
  74. PERM    *pi;
  75. {
  76.     int        i, j, k;
  77.  
  78.     if ( ! pi )
  79.     error(E_NULL,"px_rand");
  80.  
  81.     for ( i = 0; i < 3*pi->size; i++ )
  82.     {
  83.     j = (rand() >> 8) % pi->size;
  84.     k = (rand() >> 8) % pi->size;
  85.     px_transp(pi,j,k);
  86.     }
  87.  
  88.     return pi;
  89. }
  90. #ifdef RISC_OS
  91. #define SAVE_FILE "asx521_mat"
  92. #else
  93. #define    SAVE_FILE    "asx5213a.mat"
  94. #endif
  95. #define    MATLAB_NAME    "alpha"
  96. char    name[81] = MATLAB_NAME;
  97.  
  98. int main(argc, argv)
  99. int    argc;
  100. char    *argv[];
  101. {
  102.    VEC    *x = VNULL, *y = VNULL, *z = VNULL, *u = VNULL, *v = VNULL, 
  103.         *w = VNULL;
  104.    VEC    *diag = VNULL, *beta = VNULL;
  105.    PERM    *pi1 = PNULL, *pi2 = PNULL, *pi3 = PNULL, *pivot = PNULL, 
  106.         *blocks = PNULL;
  107.    MAT    *A = MNULL, *Atmp = MNULL, *B = MNULL, *C = MNULL, *D = MNULL,
  108.         *Q = MNULL, *Qtmp = MNULL, *U = MNULL;
  109.    BAND *bA, *bB, *bC;
  110.    Real    cond_est, s1, s2, s3;
  111.    int    i, j, seed;
  112.    FILE    *fp;
  113.    char    *cp;
  114.  
  115.  
  116.     mem_info_on(TRUE);
  117.  
  118.     setbuf(stdout,(char *)NULL);
  119.  
  120.     seed = 1111;
  121.     if ( argc > 2 )
  122.     {
  123.     printf("usage: %s [seed]\n",argv[0]);
  124.     exit(0);
  125.     }
  126.     else if ( argc == 2 )
  127.     sscanf(argv[1], "%d", &seed);
  128.  
  129.     /* set seed for rand() */
  130.     smrand(seed);
  131.  
  132.     mem_stat_mark(1);
  133.  
  134.     /* print version information */
  135.     m_version();
  136.  
  137.     printf("# grep \"^Error\" the output for a listing of errors\n");
  138.     printf("# Don't panic if you see \"Error\" appearing; \n");
  139.     printf("# Also check the reported size of error\n");
  140.     printf("# This program uses randomly generated problems and therefore\n");
  141.     printf("# may occasionally produce ill-conditioned problems\n");
  142.     printf("# Therefore check the size of the error compared with MACHEPS\n");
  143.     printf("# If the error is within 1000*MACHEPS then don't worry\n");
  144.     printf("# If you get an error of size 0.1 or larger there is \n");
  145.     printf("# probably a bug in the code or the compilation procedure\n\n");
  146.     printf("# seed = %d\n",seed);
  147.  
  148.     printf("# Check: MACHEPS = %g\n",MACHEPS);
  149.     /* allocate, initialise, copy and resize operations */
  150.     /* VEC */
  151.     notice("vector initialise, copy & resize");
  152.     x = v_get(12);
  153.     y = v_get(15);
  154.     z = v_get(12);
  155.     v_rand(x);
  156.     v_rand(y);
  157.     z = v_copy(x,z);
  158.     if ( v_norm2(v_sub(x,z,z)) >= MACHEPS )
  159.     errmesg("VEC copy");
  160.     v_copy(x,y);
  161.     x = v_resize(x,10);
  162.     y = v_resize(y,10);
  163.     if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
  164.     errmesg("VEC copy/resize");
  165.     x = v_resize(x,15);
  166.     y = v_resize(y,15);
  167.     if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
  168.     errmesg("VEC resize");
  169.  
  170.     /* MAT */
  171.     notice("matrix initialise, copy & resize");
  172.     A = m_get(8,5);
  173.     B = m_get(3,9);
  174.     C = m_get(8,5);
  175.     m_rand(A);
  176.     m_rand(B);
  177.     C = m_copy(A,C);
  178.     if ( m_norm_inf(m_sub(A,C,C)) >= MACHEPS )
  179.     errmesg("MAT copy");
  180.     m_copy(A,B);
  181.     A = m_resize(A,3,5);
  182.     B = m_resize(B,3,5);
  183.     if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
  184.     errmesg("MAT copy/resize");
  185.     A = m_resize(A,10,10);
  186.     B = m_resize(B,10,10);
  187.     if ( m_norm_inf(m_sub(A,B,C)) >= MACHEPS )
  188.     errmesg("MAT resize");
  189.  
  190.     MEMCHK();
  191.  
  192.     /* PERM */
  193.     notice("permutation initialise, inverting & permuting vectors");
  194.     pi1 = px_get(15);
  195.     pi2 = px_get(12);
  196.     px_rand(pi1);
  197.     v_rand(x);
  198.     px_vec(pi1,x,z);
  199.     y = v_resize(y,x->dim);
  200.     pxinv_vec(pi1,z,y);
  201.     if ( v_norm2(v_sub(x,y,z)) >= MACHEPS )
  202.     errmesg("PERMute vector");
  203.     pi2 = px_inv(pi1,pi2);
  204.     pi3 = px_mlt(pi1,pi2,PNULL);
  205.     for ( i = 0; i < pi3->size; i++ )
  206.     if ( pi3->pe[i] != i )
  207.         errmesg("PERM inverse/multiply");
  208.  
  209.     /* testing catch() etc */
  210.     notice("error handling routines");
  211.     catch(E_NULL,
  212.       catchall(v_add(VNULL,VNULL,VNULL);
  213.              errmesg("tracecatch() failure"),
  214.              printf("# tracecatch() caught error\n");
  215.              error(E_NULL,"main"));
  216.                  errmesg("catch() failure"),
  217.       printf("# catch() caught E_NULL error\n"));
  218.  
  219.     /* testing attaching a new error list (error list 2) */
  220.  
  221.     notice("attaching error lists");
  222.     printf("# IT IS NOT A REAL WARNING ... \n");
  223.     err_list_attach(2,MAX_TEST_ERR,test_err_list,TRUE);
  224.     if (!err_is_list_attached(2)) 
  225.        errmesg("attaching the error list 2");
  226.     ev_err(__FILE__,1,__LINE__,"main",2);
  227.     err_list_free(2);
  228.     if (err_is_list_attached(2)) 
  229.        errmesg("detaching the error list 2");
  230.  
  231.     /* testing inner products and v_mltadd() etc */
  232.     notice("inner products and linear combinations");
  233.     u = v_get(x->dim);
  234.     v_rand(u);
  235.     v_rand(x);
  236.     v_resize(y,x->dim);
  237.     v_rand(y);
  238.     v_mltadd(y,x,-in_prod(x,y)/in_prod(x,x),z);
  239.     if ( fabs(in_prod(x,z)) >= MACHEPS*x->dim )
  240.     errmesg("v_mltadd()/in_prod()");
  241.     s1 = -in_prod(x,y)/(v_norm2(x)*v_norm2(x));
  242.     sv_mlt(s1,x,u);
  243.     v_add(y,u,u);
  244.     if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
  245.     errmesg("sv_mlt()/v_norm2()");
  246.  
  247. #ifdef ANSI_C 
  248.     v_linlist(u,x,s1,y,1.0,VNULL);
  249.     if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
  250.     errmesg("v_linlist()");
  251. #endif
  252. #ifdef VARARGS
  253.     v_linlist(u,x,s1,y,1.0,VNULL);
  254.     if ( v_norm2(v_sub(u,z,u)) >= MACHEPS*x->dim )
  255.     errmesg("v_linlist()");
  256. #endif
  257.  
  258.  
  259.     MEMCHK();
  260.  
  261.     /* vector norms */
  262.     notice("vector norms");
  263.     x = v_resize(x,12);
  264.     v_rand(x);
  265.     for ( i = 0; i < x->dim; i++ )
  266.     if ( v_entry(x,i) >= 0.5 )
  267.         v_set_val(x,i,1.0);
  268.         else
  269.         v_set_val(x,i,-1.0);
  270.     s1 = v_norm1(x);
  271.     s2 = v_norm2(x);    
  272.     s3 = v_norm_inf(x);
  273.     if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
  274.      fabs(s2 - sqrt((Real)(x->dim))) >= MACHEPS*x->dim ||
  275.      fabs(s3 - 1.0) >= MACHEPS )
  276.     errmesg("v_norm1/2/_inf()");
  277.  
  278.     /* test matrix multiply etc */
  279.     notice("matrix multiply and invert");
  280.     A = m_resize(A,10,10);
  281.     B = m_resize(B,10,10);
  282.     m_rand(A);
  283.     m_inverse(A,B);
  284.     m_mlt(A,B,C);
  285.     for ( i = 0; i < C->m; i++ )
  286.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  287.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  288.     errmesg("m_inverse()/m_mlt()");
  289.  
  290.     MEMCHK();
  291.  
  292.     /* ... and transposes */
  293.     notice("transposes and transpose-multiplies");
  294.     m_transp(A,A);    /* can do square matrices in situ */
  295.     mtrm_mlt(A,B,C);
  296.     for ( i = 0; i < C->m; i++ )
  297.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  298.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  299.     errmesg("m_transp()/mtrm_mlt()");
  300.     m_transp(A,A);
  301.     m_transp(B,B);
  302.     mmtr_mlt(A,B,C);
  303.     for ( i = 0; i < C->m; i++ )
  304.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  305.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  306.     errmesg("m_transp()/mmtr_mlt()");
  307.     sm_mlt(3.71,B,B);
  308.     mmtr_mlt(A,B,C);
  309.     for ( i = 0; i < C->m; i++ )
  310.     m_set_val(C,i,i,m_entry(C,i,i)-3.71);
  311.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  312.     errmesg("sm_mlt()/mmtr_mlt()");
  313.     m_transp(B,B);
  314.     sm_mlt(1.0/3.71,B,B);
  315.  
  316.     MEMCHK();
  317.  
  318.     /* ... and matrix-vector multiplies */
  319.     notice("matrix-vector multiplies");
  320.     x = v_resize(x,A->n);
  321.     y = v_resize(y,A->m);
  322.     z = v_resize(z,A->m);
  323.     u = v_resize(u,A->n);
  324.     v_rand(x);
  325.     v_rand(y);
  326.     mv_mlt(A,x,z);
  327.     s1 = in_prod(y,z);
  328.     vm_mlt(A,y,u);
  329.     s2 = in_prod(u,x);
  330.     if ( fabs(s1 - s2) >= (MACHEPS*x->dim)*x->dim )
  331.     errmesg("mv_mlt()/vm_mlt()");
  332.     mv_mlt(B,z,u);
  333.     if ( v_norm2(v_sub(u,x,u)) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  334.     errmesg("mv_mlt()/m_inverse()");
  335.  
  336.     MEMCHK();
  337.  
  338.     /* get/set row/col */
  339.     notice("getting and setting rows and cols");
  340.     x = v_resize(x,A->n);
  341.     y = v_resize(y,B->m);
  342.     x = get_row(A,3,x);
  343.     y = get_col(B,3,y);
  344.     if ( fabs(in_prod(x,y) - 1.0) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  345.     errmesg("get_row()/get_col()");
  346.     sv_mlt(-1.0,x,x);
  347.     sv_mlt(-1.0,y,y);
  348.     set_row(A,3,x);
  349.     set_col(B,3,y);
  350.     m_mlt(A,B,C);
  351.     for ( i = 0; i < C->m; i++ )
  352.     m_set_val(C,i,i,m_entry(C,i,i)-1.0);
  353.     if ( m_norm_inf(C) >= MACHEPS*m_norm_inf(A)*m_norm_inf(B)*5 )
  354.     errmesg("set_row()/set_col()");
  355.  
  356.     MEMCHK();
  357.  
  358.     /* matrix norms */
  359.     notice("matrix norms");
  360.     A = m_resize(A,11,15);
  361.     m_rand(A);
  362.     s1 = m_norm_inf(A);
  363.     B = m_transp(A,B);
  364.     s2 = m_norm1(B);
  365.     if ( fabs(s1 - s2) >= MACHEPS*A->m )
  366.     errmesg("m_norm1()/m_norm_inf()");
  367.     C = mtrm_mlt(A,A,C);
  368.     s1 = 0.0;
  369.     for ( i = 0; i < C->m && i < C->n; i++ )
  370.     s1 += m_entry(C,i,i);
  371.     if ( fabs(sqrt(s1) - m_norm_frob(A)) >= MACHEPS*A->m*A->n )
  372.     errmesg("m_norm_frob");
  373.  
  374.     MEMCHK();
  375.     
  376.     /* permuting rows and columns */
  377.     notice("permuting rows & cols");
  378.     A = m_resize(A,11,15);
  379.     B = m_resize(B,11,15);
  380.     pi1 = px_resize(pi1,A->m);
  381.     px_rand(pi1);
  382.     x = v_resize(x,A->n);
  383.     y = mv_mlt(A,x,y);
  384.     px_rows(pi1,A,B);
  385.     px_vec(pi1,y,z);
  386.     mv_mlt(B,x,u);
  387.     if ( v_norm2(v_sub(z,u,u)) >= MACHEPS*A->m )
  388.     errmesg("px_rows()");
  389.     pi1 = px_resize(pi1,A->n);
  390.     px_rand(pi1);
  391.     px_cols(pi1,A,B);
  392.     pxinv_vec(pi1,x,z);
  393.     mv_mlt(B,z,u);
  394.     if ( v_norm2(v_sub(y,u,u)) >= MACHEPS*A->n )
  395.     errmesg("px_cols()");
  396.  
  397.     MEMCHK();
  398.  
  399.     /* MATLAB save/load */
  400.     notice("MATLAB save/load");
  401.     A = m_resize(A,12,11);
  402.     if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
  403.     printf("Cannot perform MATLAB save/load test\n");
  404.     else
  405.     {
  406.     m_rand(A);
  407.     m_save(fp, A, name);
  408.     fclose(fp);
  409.     if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
  410.         printf("Cannot open save file \"%s\"\n",SAVE_FILE);
  411.     else
  412.     {
  413.         M_FREE(B);
  414.         B = m_load(fp,&cp);
  415.         if ( strcmp(name,cp) || m_norm1(m_sub(A,B,B)) >= MACHEPS*A->m )
  416.         errmesg("mload()/m_save()");
  417.     }
  418.     }
  419.  
  420.     MEMCHK();
  421.  
  422.     /* Now, onto matrix factorisations */
  423.     A = m_resize(A,10,10);
  424.     B = m_resize(B,A->m,A->n);
  425.     m_copy(A,B);
  426.     x = v_resize(x,A->n);
  427.     y = v_resize(y,A->m);
  428.     z = v_resize(z,A->n);
  429.     u = v_resize(u,A->m);
  430.     v_rand(x);
  431.     mv_mlt(B,x,y);
  432.     z = v_copy(x,z);
  433.  
  434.     notice("LU factor/solve");
  435.     pivot = px_get(A->m);
  436.     LUfactor(A,pivot);
  437.     tracecatch(LUsolve(A,pivot,y,x),"main");
  438.     tracecatch(cond_est = LUcondest(A,pivot),"main");
  439.     printf("# cond(A) approx= %g\n", cond_est);
  440.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
  441.     {
  442.     errmesg("LUfactor()/LUsolve()");
  443.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  444.            v_norm2(v_sub(x,z,u)), MACHEPS);
  445.     }
  446.  
  447.     v_copy(y,x);
  448.     tracecatch(LUsolve(A,pivot,x,x),"main");
  449.     tracecatch(cond_est = LUcondest(A,pivot),"main");
  450.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
  451.     {
  452.     errmesg("LUfactor()/LUsolve()");
  453.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  454.            v_norm2(v_sub(x,z,u)), MACHEPS);
  455.     }
  456.  
  457.     vm_mlt(B,z,y);
  458.     v_copy(y,x);
  459.     tracecatch(LUTsolve(A,pivot,x,x),"main");
  460.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est)
  461.     {
  462.     errmesg("LUfactor()/LUTsolve()");
  463.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  464.            v_norm2(v_sub(x,z,u)), MACHEPS);
  465.     }
  466.  
  467.     MEMCHK();
  468.  
  469.     /* QR factorisation */
  470.     m_copy(B,A);
  471.     mv_mlt(B,z,y);
  472.     notice("QR factor/solve:");
  473.     diag = v_get(A->m);
  474.     beta = v_get(A->m);
  475.     QRfactor(A,diag);
  476.     QRsolve(A,diag,y,x);
  477.     if ( v_norm2(v_sub(x,z,u)) >= MACHEPS*v_norm2(x)*cond_est )
  478.     {
  479.     errmesg("QRfactor()/QRsolve()");
  480.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  481.            v_norm2(v_sub(x,z,u)), MACHEPS);
  482.     }
  483.     Q = m_get(A->m,A->m);
  484.     makeQ(A,diag,Q);
  485.     makeR(A,A);
  486.     m_mlt(Q,A,C);
  487.     m_sub(B,C,C);
  488.     if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  489.     {
  490.     errmesg("QRfactor()/makeQ()/makeR()");
  491.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  492.            m_norm1(C), MACHEPS);
  493.     }
  494.  
  495.     MEMCHK();
  496.  
  497.     /* now try with a non-square matrix */
  498.     A = m_resize(A,15,7);
  499.     m_rand(A);
  500.     B = m_copy(A,B);
  501.     diag = v_resize(diag,A->n);
  502.     beta = v_resize(beta,A->n);
  503.     x = v_resize(x,A->n);
  504.     y = v_resize(y,A->m);
  505.     v_rand(y);
  506.     QRfactor(A,diag);
  507.     x = QRsolve(A,diag,y,x);
  508.     /* z is the residual vector */
  509.     mv_mlt(B,x,z);    v_sub(z,y,z);
  510.     /* check B^T.z = 0 */
  511.     vm_mlt(B,z,u);
  512.     if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
  513.     {
  514.     errmesg("QRfactor()/QRsolve()");
  515.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  516.            v_norm2(u), MACHEPS);
  517.     }
  518.     Q = m_resize(Q,A->m,A->m);
  519.     makeQ(A,diag,Q);
  520.     makeR(A,A);
  521.     m_mlt(Q,A,C);
  522.     m_sub(B,C,C);
  523.     if ( m_norm1(C) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  524.     {
  525.     errmesg("QRfactor()/makeQ()/makeR()");
  526.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  527.            m_norm1(C), MACHEPS);
  528.     }
  529.     D = m_get(A->m,Q->m);
  530.     mtrm_mlt(Q,Q,D);
  531.     for ( i = 0; i < D->m; i++ )
  532.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  533.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
  534.     {
  535.     errmesg("QRfactor()/makeQ()/makeR()");
  536.     printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
  537.            m_norm1(D), MACHEPS);
  538.     }
  539.  
  540.     MEMCHK();
  541.  
  542.     /* QRCP factorisation */
  543.     m_copy(B,A);
  544.     notice("QR factor/solve with column pivoting");
  545.     pivot = px_resize(pivot,A->n);
  546.     QRCPfactor(A,diag,pivot);
  547.     z = v_resize(z,A->n);
  548.     QRCPsolve(A,diag,pivot,y,z);
  549.     /* pxinv_vec(pivot,z,x); */
  550.     /* now compute residual (z) vector */
  551.     mv_mlt(B,x,z);    v_sub(z,y,z);
  552.     /* check B^T.z = 0 */
  553.     vm_mlt(B,z,u);
  554.     if ( v_norm2(u) >= MACHEPS*m_norm1(B)*v_norm2(y) )
  555.     {
  556.     errmesg("QRCPfactor()/QRsolve()");
  557.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  558.            v_norm2(u), MACHEPS);
  559.     }
  560.  
  561.     Q = m_resize(Q,A->m,A->m);
  562.     makeQ(A,diag,Q);
  563.     makeR(A,A);
  564.     m_mlt(Q,A,C);
  565.     M_FREE(D);
  566.     D = m_get(B->m,B->n);
  567.     px_cols(pivot,C,D);
  568.     m_sub(B,D,D);
  569.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm1(B) )
  570.     {
  571.     errmesg("QRCPfactor()/makeQ()/makeR()");
  572.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  573.            m_norm1(D), MACHEPS);
  574.     }
  575.  
  576.     MEMCHK();
  577.  
  578.     /* Cholesky and LDL^T factorisation */
  579.     /* Use these for normal equations approach */
  580.     notice("Cholesky factor/solve");
  581.     mtrm_mlt(B,B,A);
  582.     CHfactor(A);
  583.     u = v_resize(u,B->n);
  584.     vm_mlt(B,y,u);
  585.     z = v_resize(z,B->n);
  586.     CHsolve(A,u,z);
  587.     v_sub(x,z,z);
  588.     if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
  589.     {
  590.     errmesg("CHfactor()/CHsolve()");
  591.     printf("# Cholesky solution error = %g [cf MACHEPS = %g]\n",
  592.            v_norm2(z), MACHEPS);
  593.     }
  594.     /* modified Cholesky factorisation should be identical with Cholesky
  595.        factorisation provided the matrix is "sufficiently positive definite */
  596.     mtrm_mlt(B,B,C);
  597.     MCHfactor(C,MACHEPS);
  598.     m_sub(A,C,C);
  599.     if ( m_norm1(C) >= MACHEPS*m_norm1(A) )
  600.     {
  601.     errmesg("MCHfactor()");
  602.     printf("# Modified Cholesky error = %g [cf MACHEPS = %g]\n",
  603.            m_norm1(C), MACHEPS);
  604.     }
  605.     /* now test the LDL^T factorisation -- using a negative def. matrix */
  606.     mtrm_mlt(B,B,A);
  607.     sm_mlt(-1.0,A,A);
  608.     m_copy(A,C);
  609.     LDLfactor(A);
  610.     LDLsolve(A,u,z);
  611.     w = v_get(A->m);
  612.     mv_mlt(C,z,w);
  613.     v_sub(w,u,w);
  614.     if ( v_norm2(w) >= MACHEPS*v_norm2(u)*m_norm1(C) )
  615.     {
  616.     errmesg("LDLfactor()/LDLsolve()");
  617.     printf("# LDL^T residual = %g [cf MACHEPS = %g]\n",
  618.            v_norm2(w), MACHEPS);
  619.     }
  620.     v_add(x,z,z);
  621.     if ( v_norm2(z) >= MACHEPS*v_norm2(x)*100 )
  622.     {
  623.     errmesg("LDLfactor()/LDLsolve()");
  624.     printf("# LDL^T solution error = %g [cf MACHEPS = %g]\n",
  625.            v_norm2(z), MACHEPS);
  626.     }
  627.  
  628.     MEMCHK();
  629.  
  630.     /* and now the Bunch-Kaufman-Parlett method */
  631.     /* set up D to be an indefinite diagonal matrix */
  632.     notice("Bunch-Kaufman-Parlett factor/solve");
  633.  
  634.     D = m_resize(D,B->m,B->m);
  635.     m_zero(D);
  636.     w = v_resize(w,B->m);
  637.     v_rand(w);
  638.     for ( i = 0; i < w->dim; i++ )
  639.     if ( v_entry(w,i) >= 0.5 )
  640.         m_set_val(D,i,i,1.0);
  641.     else
  642.         m_set_val(D,i,i,-1.0);
  643.     /* set A <- B^T.D.B */
  644.     C = m_resize(C,B->n,B->n);
  645.     C = mtrm_mlt(B,D,C);
  646.     A = m_mlt(C,B,A);
  647.     C = m_resize(C,B->n,B->n);
  648.     C = m_copy(A,C);
  649.     /* ... and use BKPfactor() */
  650.     blocks = px_get(A->m);
  651.     pivot = px_resize(pivot,A->m);
  652.     x = v_resize(x,A->m);
  653.     y = v_resize(y,A->m);
  654.     z = v_resize(z,A->m);
  655.     v_rand(x);
  656.     mv_mlt(A,x,y);
  657.     BKPfactor(A,pivot,blocks);
  658.     printf("# BKP pivot =\n");    px_output(pivot);
  659.     printf("# BKP blocks =\n");    px_output(blocks);
  660.     BKPsolve(A,pivot,blocks,y,z);
  661.     /* compute & check residual */
  662.     mv_mlt(C,z,w);
  663.     v_sub(w,y,w);
  664.     if ( v_norm2(w) >= MACHEPS*m_norm1(C)*v_norm2(z) )
  665.     {
  666.     errmesg("BKPfactor()/BKPsolve()");
  667.     printf("# BKP residual size = %g [cf MACHEPS = %g]\n",
  668.            v_norm2(w), MACHEPS);
  669.     }
  670.  
  671.     /* check update routines */
  672.     /* check LDLupdate() first */
  673.     notice("update L.D.L^T routine");
  674.     A = mtrm_mlt(B,B,A);
  675.     m_resize(C,A->m,A->n);
  676.     C = m_copy(A,C);
  677.     LDLfactor(A);
  678.     s1 = 3.7;
  679.     w = v_resize(w,A->m);
  680.     v_rand(w);
  681.     for ( i = 0; i < C->m; i++ )
  682.     for ( j = 0; j < C->n; j++ )
  683.         m_set_val(C,i,j,m_entry(C,i,j)+s1*v_entry(w,i)*v_entry(w,j));
  684.     LDLfactor(C);
  685.     LDLupdate(A,w,s1);
  686.     /* zero out strictly upper triangular parts of A and C */
  687.     for ( i = 0; i < A->m; i++ )
  688.     for ( j = i+1; j < A->n; j++ )
  689.     {
  690.         m_set_val(A,i,j,0.0);
  691.         m_set_val(C,i,j,0.0);
  692.     }
  693.     if ( m_norm1(m_sub(A,C,C)) >= sqrt(MACHEPS)*m_norm1(A) )
  694.     {
  695.     errmesg("LDLupdate()");
  696.     printf("# LDL update matrix error = %g [cf MACHEPS = %g]\n",
  697.            m_norm1(C), MACHEPS);
  698.     }
  699.  
  700.  
  701.     /* BAND MATRICES */
  702.  
  703. #define COL 40
  704. #define UDIAG  5
  705. #define LDIAG  2
  706.  
  707.    smrand(101);
  708.    bA = bd_get(LDIAG,UDIAG,COL);
  709.    bB = bd_get(LDIAG,UDIAG,COL);
  710.    bC = bd_get(LDIAG,UDIAG,COL);
  711.    A = m_resize(A,COL,COL);
  712.    B = m_resize(B,COL,COL);
  713.    pivot = px_resize(pivot,COL);
  714.    x = v_resize(x,COL);
  715.    w = v_resize(w,COL);
  716.    z = v_resize(z,COL);
  717.  
  718.    m_rand(A); 
  719.    /* generate band matrix */
  720.    mat2band(A,LDIAG,UDIAG,bA);
  721.    band2mat(bA,A);    /* now A is banded */
  722.    bB = bd_copy(bA,bB); 
  723.  
  724.    v_rand(x);  
  725.    mv_mlt(A,x,w);
  726.    z = v_copy(w,z);
  727.  
  728.    notice("band LU factorization");
  729.    bdLUfactor(bA,pivot);
  730.  
  731.    /* pivot will be changed */
  732.    bdLUsolve(bA,pivot,z,z);
  733.    v_sub(x,z,z);
  734.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  735.       errmesg("incorrect solution (band LU factorization)");
  736.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  737.          v_norm2(z),MACHEPS);
  738.    }
  739.  
  740.    /* solve transpose system */
  741.  
  742.    notice("band LU factorization for transpose system");
  743.    m_transp(A,B);
  744.    mv_mlt(B,x,w);
  745.  
  746.    bd_copy(bB,bA);
  747.    bd_transp(bA,bA);  
  748.    /* transposition in situ */
  749.    bd_transp(bA,bA);
  750.    bd_transp(bA,bB);
  751.  
  752.    bdLUfactor(bB,pivot);
  753.  
  754.    bdLUsolve(bB,pivot,w,z);
  755.    v_sub(x,z,z);
  756.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  757.       errmesg("incorrect solution (band transposed LU factorization)");
  758.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  759.          v_norm2(z),MACHEPS);
  760.    }
  761.  
  762.  
  763.    /* Cholesky factorization */
  764.  
  765.    notice("band Choleski LDL' factorization");
  766.    m_add(A,B,A);  /* symmetric matrix */
  767.    for (i=0; i < COL; i++)     /* positive definite */
  768.      A->me[i][i] += 2*LDIAG;   
  769.  
  770.    mat2band(A,LDIAG,LDIAG,bA);
  771.    band2mat(bA,A);              /* corresponding matrix A */
  772.  
  773.    v_rand(x);
  774.    mv_mlt(A,x,w);
  775.    z = v_copy(w,z);
  776.    
  777.    bdLDLfactor(bA);
  778.  
  779.    z = bdLDLsolve(bA,z,z);
  780.    v_sub(x,z,z);
  781.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  782.       errmesg("incorrect solution (band LDL' factorization)");
  783.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  784.          v_norm2(z),MACHEPS);
  785.    }
  786.  
  787.    /* new bandwidths */
  788.    m_rand(A);
  789.    bA = bd_resize(bA,UDIAG,LDIAG,COL);
  790.    bB = bd_resize(bB,UDIAG,LDIAG,COL);
  791.    mat2band(A,UDIAG,LDIAG,bA);
  792.    band2mat(bA,A);
  793.    bd_copy(bA,bB);
  794.  
  795.    mv_mlt(A,x,w);
  796.  
  797.    notice("band LU factorization (resized)");
  798.    bdLUfactor(bA,pivot);
  799.  
  800.    /* pivot will be changed */
  801.    bdLUsolve(bA,pivot,w,z);
  802.    v_sub(x,z,z);
  803.    if (v_norm2(z) > v_norm2(x)*sqrt(MACHEPS)) {
  804.       errmesg("incorrect solution (band LU factorization)");
  805.       printf(" ||exact sol. - computed sol.|| = %g [MACHEPS = %g]\n",
  806.          v_norm2(z),MACHEPS);
  807.    }
  808.  
  809.    /* testing transposition */
  810.  
  811.    notice("band matrix transposition");
  812.    m_zero(bA->mat);
  813.    bd_copy(bB,bA);
  814.    m_zero(bB->mat);
  815.    bd_copy(bA,bB);
  816.  
  817.    bd_transp(bB,bB);
  818.    bd_transp(bB,bB);
  819.  
  820.    m_zero(bC->mat);
  821.    bd_copy(bB,bC);
  822.  
  823.    m_sub(bA->mat,bC->mat,bC->mat);
  824.    if (m_norm_inf(bC->mat) > MACHEPS*bC->mat->n) {
  825.       errmesg("band transposition");
  826.       printf(" difference ||A - (A')'|| = %g\n",m_norm_inf(bC->mat));
  827.    }
  828.  
  829.    bd_free(bA);
  830.    bd_free(bB);
  831.    bd_free(bC);
  832.  
  833.  
  834.     MEMCHK();
  835.  
  836.     /* now check QRupdate() routine */
  837.     notice("update QR routine");
  838.  
  839.     B = m_resize(B,15,7);
  840.     A = m_resize(A,B->m,B->n);
  841.     m_copy(B,A);
  842.     diag = v_resize(diag,A->n);
  843.     beta = v_resize(beta,A->n);
  844.     QRfactor(A,diag);
  845.     Q = m_resize(Q,A->m,A->m);
  846.     makeQ(A,diag,Q);
  847.     makeR(A,A);
  848.     m_resize(C,A->m,A->n);
  849.     w = v_resize(w,A->m);
  850.     v = v_resize(v,A->n);
  851.     u = v_resize(u,A->m);
  852.     v_rand(w);
  853.     v_rand(v);
  854.     vm_mlt(Q,w,u);
  855.     QRupdate(Q,A,u,v);
  856.     m_mlt(Q,A,C);
  857.     for ( i = 0; i < B->m; i++ )
  858.     for ( j = 0; j < B->n; j++ )
  859.         m_set_val(B,i,j,m_entry(B,i,j)+v_entry(w,i)*v_entry(v,j));
  860.     m_sub(B,C,C);
  861.     if ( m_norm1(C) >= MACHEPS*m_norm1(A)*m_norm1(Q)*2 )
  862.     {
  863.     errmesg("QRupdate()");
  864.     printf("# Reconstruction error in QR update = %g [cf MACHEPS = %g]\n",
  865.            m_norm1(C), MACHEPS);
  866.     }
  867.     m_resize(D,Q->m,Q->n);
  868.     mtrm_mlt(Q,Q,D);
  869.     for ( i = 0; i < D->m; i++ )
  870.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  871.     if ( m_norm1(D) >= 10*MACHEPS*m_norm1(Q)*m_norm_inf(Q) )
  872.     {
  873.     errmesg("QRupdate()");
  874.     printf("# QR update orthogonality error = %g [cf MACHEPS = %g]\n",
  875.            m_norm1(D), MACHEPS);
  876.     }
  877.  
  878.     /* Now check eigenvalue/SVD routines */
  879.     notice("eigenvalue and SVD routines");
  880.     A = m_resize(A,11,11);
  881.     B = m_resize(B,A->m,A->n);
  882.     C = m_resize(C,A->m,A->n);
  883.     D = m_resize(D,A->m,A->n);
  884.     Q = m_resize(Q,A->m,A->n);
  885.  
  886.     m_rand(A);
  887.     /* A <- A + A^T  for symmetric case */
  888.     m_add(A,m_transp(A,C),A);
  889.     u = v_resize(u,A->m);
  890.     u = symmeig(A,Q,u);
  891.     m_zero(B);
  892.     for ( i = 0; i < B->m; i++ )
  893.     m_set_val(B,i,i,v_entry(u,i));
  894.     m_mlt(Q,B,C);
  895.     mmtr_mlt(C,Q,D);
  896.     m_sub(A,D,D);
  897.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*v_norm_inf(u)*3 )
  898.     {
  899.     errmesg("symmeig()");
  900.     printf("# Reconstruction error = %g [cf MACHEPS = %g]\n",
  901.            m_norm1(D), MACHEPS);
  902.     }
  903.     mtrm_mlt(Q,Q,D);
  904.     for ( i = 0; i < D->m; i++ )
  905.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  906.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*3 )
  907.     {
  908.     errmesg("symmeig()");
  909.     printf("# symmeig() orthogonality error = %g [cf MACHEPS = %g]\n",
  910.            m_norm1(D), MACHEPS);
  911.     }
  912.  
  913.     MEMCHK();
  914.  
  915.     /* now test (real) Schur decomposition */
  916.     /* m_copy(A,B); */
  917.     M_FREE(A);
  918.     A = m_get(11,11);
  919.     m_rand(A);
  920.     B = m_copy(A,B);
  921.     MEMCHK();
  922.  
  923.     B = schur(B,Q);
  924.     MEMCHK();
  925.  
  926.     m_mlt(Q,B,C);
  927.     mmtr_mlt(C,Q,D);
  928.     MEMCHK();
  929.     m_sub(A,D,D);
  930.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*m_norm1(B)*5 )
  931.     {
  932.     errmesg("schur()");
  933.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  934.            m_norm1(D), MACHEPS);
  935.     }
  936.  
  937.     /* orthogonality check */
  938.     mmtr_mlt(Q,Q,D);
  939.     for ( i = 0; i < D->m; i++ )
  940.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  941.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*10 )
  942.     {
  943.     errmesg("schur()");
  944.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  945.            m_norm1(D), MACHEPS);
  946.     }
  947.  
  948.     MEMCHK();
  949.  
  950.     /* check for termination */
  951.     Atmp = m_get(2,2);
  952.     Qtmp = m_get(2,2);
  953.     /* this is a 2 x 2 matrix with real eigenvalues */
  954.     Atmp->me[0][0] = 1;
  955.     Atmp->me[1][1] = 1;
  956.     Atmp->me[0][1] = 4;
  957.     Atmp->me[1][0] = 1;
  958.     schur(Atmp,Qtmp);
  959.  
  960.     MEMCHK();
  961.  
  962.     /* now test SVD */
  963.     A = m_resize(A,11,7);
  964.     m_rand(A);
  965.     U = m_get(A->n,A->n);
  966.     Q = m_resize(Q,A->m,A->m);
  967.     u = v_resize(u,max(A->m,A->n));
  968.     svd(A,Q,U,u);
  969.     /* check reconstruction of A */
  970.     D = m_resize(D,A->m,A->n);
  971.     C = m_resize(C,A->m,A->n);
  972.     m_zero(D);
  973.     for ( i = 0; i < min(A->m,A->n); i++ )
  974.     m_set_val(D,i,i,v_entry(u,i));
  975.     mtrm_mlt(Q,D,C);
  976.     m_mlt(C,U,D);
  977.     m_sub(A,D,D);
  978.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(Q)*m_norm1(A) )
  979.     {
  980.     errmesg("svd()");
  981.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  982.            m_norm1(D), MACHEPS);
  983.     }
  984.     /* check orthogonality of Q and U */
  985.     D = m_resize(D,Q->n,Q->n);
  986.     mtrm_mlt(Q,Q,D);
  987.     for ( i = 0; i < D->m; i++ )
  988.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  989.     if ( m_norm1(D) >= MACHEPS*m_norm1(Q)*m_norm_inf(Q)*5 )
  990.     {
  991.     errmesg("svd()");
  992.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  993.            m_norm1(D), MACHEPS);
  994.     }
  995.     D = m_resize(D,U->n,U->n);
  996.     mtrm_mlt(U,U,D);
  997.     for ( i = 0; i < D->m; i++ )
  998.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  999.     if ( m_norm1(D) >= MACHEPS*m_norm1(U)*m_norm_inf(U)*5 )
  1000.     {
  1001.     errmesg("svd()");
  1002.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  1003.            m_norm1(D), MACHEPS);
  1004.     }
  1005.     for ( i = 0; i < u->dim; i++ )
  1006.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  1007.                   v_entry(u,i+1) > v_entry(u,i)) )
  1008.         break;
  1009.     if ( i < u->dim )
  1010.     {
  1011.     errmesg("svd()");
  1012.     printf("# SVD sorting error\n");
  1013.     }
  1014.  
  1015.  
  1016.     /* test of long vectors */
  1017.     notice("Long vectors");
  1018.     x = v_resize(x,100000);
  1019.     y = v_resize(y,100000);
  1020.     z = v_resize(z,100000);
  1021.     v_rand(x);
  1022.     v_rand(y);
  1023.     v_mltadd(x,y,3.0,z);
  1024.     sv_mlt(1.0/3.0,z,z);
  1025.     v_mltadd(z,x,-1.0/3.0,z);
  1026.     v_sub(z,y,x);
  1027.     if (v_norm2(x) >= MACHEPS*(x->dim)) {
  1028.        errmesg("long vectors");
  1029.        printf(" norm = %g\n",v_norm2(x));
  1030.     }
  1031.  
  1032.     mem_stat_free(1);
  1033.  
  1034.     MEMCHK();
  1035.  
  1036.     /**************************************************
  1037.     VEC        *x, *y, *z, *u, *v, *w;
  1038.     VEC        *diag, *beta;
  1039.     PERM    *pi1, *pi2, *pi3, *pivot, *blocks;
  1040.     MAT        *A, *B, *C, *D, *Q, *U;
  1041.     **************************************************/
  1042.     V_FREE(x);        V_FREE(y);    V_FREE(z);
  1043.     V_FREE(u);        V_FREE(v);    V_FREE(w);
  1044.     V_FREE(diag);    V_FREE(beta);
  1045.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pi3);
  1046.     PX_FREE(pivot);    PX_FREE(blocks);
  1047.     M_FREE(A);        M_FREE(B);    M_FREE(C);
  1048.     M_FREE(D);        M_FREE(Q);    M_FREE(U);
  1049.     M_FREE(Atmp);    M_FREE(Qtmp);
  1050.  
  1051.     MEMCHK();
  1052.     printf("# Finished torture test\n");
  1053.     mem_info();
  1054.  
  1055.     return 0;
  1056. }
  1057.  
  1058.